FAIRE DES CARTOGRAMS DANS R
Lorem Ipsum is simply dummy text of the printing and typesetting industry. Lorem Ipsum has been the industry’s standard dummy text ever since the 1500s, when an unknown printer took a galley of type and scrambled it to make a type specimen book. It has survived not only five centuries, but also the leap into electronic typesetting, remaining essentially unchanged. It was popularised in the 1960s with the release of Letraset sheets containing Lorem Ipsum passages, and more recently with desktop publishing software like Aldus PageMaker including versions of Lorem Ipsum.
Packages
installer les packages suivants ….. + version min
install.packages(knitr)
install.packages(sf)
install.packages(mapsf)
install.packages(packcircles)
install.packages(cartogram)
install.packages(recmap)
install.packages("https://cran.r-project.org/src/contrib/Archive/cartogramR/cartogramR_1.0-1.tar.gz", repos = NULL, type = "source")Import et mise en forme des données
données INSEE (pop + CSP en 2018)
Le package sf bla bla bla…..
library(sf)
communes <- st_read("data/isere.geojson", quiet = TRUE ) %>% st_transform(2154)
data <- read.csv("data/popisrere.csv", dec = ",")
communes = merge(x = communes[,c("id","name","geometry")],
y = data[,c("id", "pop2018","agri", "art", "cadr", "interm", "emp","ouvr","retr")],
by = "id")
isere = st_union(communes)knitr::kable(communes[c(0:10),], row.names = F, digits = 1)| id | name | pop2018 | agri | art | cadr | interm | emp | ouvr | retr | geometry |
|---|---|---|---|---|---|---|---|---|---|---|
| 38001 | Les Abrets en Dauphiné | 6336 | 30.1 | 210.9 | 256.0 | 712.9 | 833.3 | 948.8 | 1397.2 | MULTIPOLYGON (((900716.2 64… |
| 38002 | Les Adrets | 1025 | 0.0 | 53.0 | 168.7 | 163.9 | 82.0 | 67.5 | 159.1 | MULTIPOLYGON (((931354.5 64… |
| 38003 | Agnin | 1127 | 0.0 | 28.9 | 72.2 | 168.6 | 139.7 | 134.9 | 192.6 | MULTIPOLYGON (((844459.1 64… |
| 38004 | L’Albenc | 1279 | 25.0 | 30.0 | 65.0 | 215.0 | 130.0 | 175.0 | 210.0 | MULTIPOLYGON (((893395.8 64… |
| 38005 | Allemond | 954 | 5.0 | 75.0 | 30.0 | 160.0 | 150.0 | 120.0 | 170.0 | MULTIPOLYGON (((934142.3 64… |
| 38006 | Allevard | 4062 | 0.0 | 176.7 | 323.2 | 535.2 | 469.6 | 469.6 | 984.2 | MULTIPOLYGON (((948125.2 64… |
| 38008 | Ambel | 28 | 10.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 10.0 | MULTIPOLYGON (((930843.9 64… |
| 38009 | Anjou | 1009 | 10.0 | 54.9 | 60.0 | 124.9 | 100.0 | 115.0 | 278.5 | MULTIPOLYGON (((847739.9 64… |
| 38010 | Annoisin-Chatelans | 686 | 14.8 | 44.4 | 54.2 | 113.3 | 64.1 | 78.9 | 113.3 | MULTIPOLYGON (((875494.4 65… |
| 38011 | Anthon | 1073 | 5.0 | 55.0 | 75.0 | 185.0 | 125.0 | 130.0 | 220.0 | MULTIPOLYGON (((866538.7 65… |
mapsf
Le package mapsf permet de faire des cartes thématiques dans R. C’est le package qui succède au package cartography.
Création d’un template cartographique
library(mapsf)
col = "#c291bc"
credits = paste0("Bronner Anne-Christine & Nicolas Lambert, 2021\n",
"Source: IGN & INSEE, 2021")
theme = mf_theme(x = "default", bg = "#f0f0f0", tab = FALSE,
pos = "center", line = 2, inner = FALSE,
fg = col, mar = c(0,0, 2, 0),cex = 1.9)
template = function(title, file, basemap = TRUE, scale = TRUE){
mf_export(
communes,
export = "png",
width = 1000,
filename = file,
res = 96,
theme = theme,
expandBB = c(-.02,0,-.02,0)
)
if (basemap == TRUE){
mf_shadow(x = communes, col = "grey50", cex = 1, add = TRUE)
mf_map(communes, col ="#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
}
mf_title(title)
if (scale == TRUE){
mf_scale(size = 20, pos = "bottomright", lwd = 1.2, cex = 1, col = "#383838", unit = "km")
}
mf_credits(
txt = credits,
pos = "bottomleft",
col = "#1a2640",
cex = 0.8,
font = 3,
bg = "#ffffff30"
)
}template("Template cartographique", "maps/template.png")
mf_map(communes, col = col, border = "white", lwd = 0.5, add = TRUE)
dev.off()Symboles proportionnels
template("Symboles proportionnels (mapsf)", "maps/prop.png")
mf_map(communes, var = "pop2018", col = col, border = "#6b4266", type = "prop",
inches = 0.8, leg_title_cex = 1.2, leg_val_cex = 0.8,
leg_title = "Nombre d'habitants, 2018")
dev.off()template("Symboles proportionnels (mapsf)", "maps/prop2.png")
mf_map(communes, var = "pop2018", col = col, border = "#6b4266", type = "prop",
inches = 0.8, leg_title_cex = 1.2, leg_val_cex = 0.8, symbol = "square",
leg_title = "Nombre d'habitants, 2018")
dev.off()Dot density maps
dotdensitymap <- function(x, var, onedot = 1, radius = 1){
x <- x[,c("id",var,"geometry")]
x[,"v"] <- round(x[,var] %>% st_drop_geometry() /onedot,0)
dots <- st_sample(x, x$v, type = "random", exact = TRUE)
circles <- st_buffer(dots, dist = radius)
return (circles)
}onedot = 500
dots = dotdensitymap(x = communes, var = "pop2018", onedot = onedot, radius = 300)
template(paste0("Un point = ",onedot," habitants"), "maps/dotdensity.png")
mf_map(dots, col = col, border = "#520a2c", lwd = 0.5, add = TRUE)
dev.off()packcircles
Le package packcirles propose 3 algorithmes simples pour déplacer des diques sur un plan 2D de telle sorte qu’ils ne se supperposent pas. Nous pouvons l’utiliser pour créer des cartogrammes de Dorling (Dorling 1996).
Création d’un ficher de données simplifié avec les coordonnées des centroides des communes.
dots = communes
st_geometry(dots) <- st_centroid(sf::st_geometry(dots),of_largest_polygon = TRUE)
dots <- data.frame(dots$id, dots["pop2018"], st_coordinates(dots))
dots = dots[,c("dots.id","X","Y","pop2018")]
colnames(dots) <- c("id","x","y","v")
dots <- dots[!is.na(dots$v),]
knitr::kable(dots[c(0:5),], row.names = F, digits = 1)| id | x | y | v |
|---|---|---|---|
| 38001 | 902026.9 | 6495943 | 6336 |
| 38002 | 934555.0 | 6466630 | 1025 |
| 38003 | 845708.9 | 6472468 | 1127 |
| 38004 | 892892.6 | 6460697 | 1279 |
| 38005 | 937029.6 | 6456880 | 954 |
La fonction circleRepelLayout() prend un ensemble de cercles dans un cadre de données et utilise la répulsion itérative pour essayer de trouver un arrangement sans chevauchement où tous les centres des cercles se trouvent à l’intérieur d’un rectangle de délimitation. Si aucun arrangement de ce type ne peut être trouvé dans le nombre maximum d’itérations spécifié, la dernière tentative est renvoyée.
library("packcircles")
k = 500 # pour ajuster la taille des cercles
itermax = 10 # nombre d'iterations
dat.init <- dots[,c("x","y","v")]
dat.init$v <- sqrt(dat.init$v * k)
simulation <- circleRepelLayout(x = dat.init, xysizecols = 1:3,
wrap = FALSE, sizetype = "radius",
maxiter = itermax, weights =1)$layout
knitr::kable(simulation[c(0:5),], row.names = F, digits = 1)| x | y | radius |
|---|---|---|
| 902120.7 | 6496140 | 1779.9 |
| 934555.0 | 6466630 | 715.9 |
| 845708.9 | 6472468 | 750.7 |
| 892892.6 | 6460697 | 799.7 |
| 937029.6 | 6456880 | 690.7 |
circles <- st_buffer(sf::st_as_sf(simulation, coords =c('x', 'y'),
crs = sf::st_crs(communes)), dist = simulation$radius)
circles$v = dots$v
template("Dorling (packcircles)", "maps/dorling1.png", scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(circles,col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()## png
## 2
cartogram
**le package cartogram est développé par Sebastian Jeworutzki. Il propose trois méthodes : Dorling, Olson et Dougenik
library(cartogram)Dorling
dorling = cartogram_dorling(communes, "pop2018", k = 1.8)template("Dorling (cartogram)", "maps/dorling2.png", scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(dorling, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()Olson
Bla bla bla cartogrammes discontinus…. (Olson 1976)
x : un objet sf weight : nom de la variable k : facteur pour augmenter la taille des polygons inplace : Si VRAI, chaque polygone est modifié à sa place initiale, si FAUX, les multi-polygones sont centrés sur leur centroïde initial
olsen <- cartogram_ncont(communes, "pop2018", k = 1, inplace = TRUE)template("Olsen (cartogram)", "maps/olsen.png", basemap = FALSE, scale = FALSE)
mf_map(isere, col = "#CCCCCC30", border = "white", lwd = 0.5, add = TRUE)
mf_map(olsen, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()Douguenik
Bla bla bla…. (Dougenik, Chrisman, and Niemeyer 1985)
x : un objet sf weight : om de la variable itermax : nombre d’itérations maximum si maxSizeError n’est pas atteint maxSizeError : la déformation s’arrete si l’erreur moyenne est inférieur à cette valeur prepare : mettre “adjust” permer d’acceler le temps de calcul. threshold : seuil pour la préparation des données
Dougenik <- cartogram_cont(communes, "pop2018", prepare = "none", itermax = 10, maxSizeError = 1.15)Calcul des érreurs
sumarea = sum(as.numeric(st_area(Dougenik)))
sumpop = sum(Dougenik$pop2018)
Dougenik$error = (as.numeric(st_area(Dougenik)) / sumarea) / (Dougenik$pop2018 / sumpop) * 100
summary(Dougenik$error)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2513 96.9756 109.4084 143.6154 130.7821 2734.5953
bks = c(min(Dougenik$error),70,80,90,100,110,120,max(Dougenik$error))
cols = c("#d53e4f", "#f46d43","#fdae61","#fee08b","#e6f598","#abdda4", "#66c2a5")Affichage de la carte
template("Dougenik (cartogram)", "maps/dougenik.png", basemap = FALSE, scale = FALSE)
mf_map(x = Dougenik, type = "choro",var = "error", pal = cols, breaks = bks, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()RecMap
Recmap (Heilmann et al. 2004) (Panse 2016) est un package développé par Christian Panse. Il perlet de concertir les géométries des unités spatiales en rectangles dont la surace est définie en fonction d’une donnée statitique sauntitative absolue. L’algorithe de RecMap est disponible ici. Développé en C++11, le package dépend des packages GA (>= 3.1), Rcpp (>= 1.0), sp(>= 1.3)
Attention, il ne fonctionne que sur un nombre limité d’unités territoriales. Et même dans ce cas, il ne fonctionne pas très bien, notamment avec les petites unités spatiales.
library(recmap)Ici un exemple sur les pays d’Europe.
europe <- st_read("data/europe.geojson")Préparation des données
coords = data.frame(st_coordinates(st_centroid(st_geometry(europe))))
bb <- lapply(st_geometry(europe), function(x){st_bbox(x)})
dx <- unlist(lapply(bb, function(x){x[3]-x[1]})) / 2
dy <- unlist(lapply(bb, function(x){x[4]-x[2]})) / 2
df <- data.frame(x = coords$X,
y = coords$Y,
dx = dx,
dy = dy,
z = europe$pop2008,
name = europe$id)knitr::kable(df[c(0:10),], row.names = F, digits = 1)| x | y | dx | dy | z | name |
|---|---|---|---|---|---|
| 4631849 | 2726828 | 281676.5 | 144612.8 | 8318592 | AT |
| 3947159 | 3073426 | 130414.4 | 108636.1 | 10666866 | BE |
| 5560592 | 2309244 | 241662.7 | 187871.0 | 7640238 | BG |
| 4190230 | 2634505 | 175597.7 | 112221.5 | 7593494 | CH |
| 6424605 | 1654151 | 89911.7 | 81559.3 | 789269 | CY |
| 4708531 | 2971173 | 247155.0 | 137212.1 | 10381130 | CZ |
| 4351269 | 3107740 | 318053.7 | 426690.0 | 82217837 | DE |
| 4321062 | 3653192 | 224577.6 | 177722.3 | 5475791 | DK |
| 5208731 | 4047869 | 171660.9 | 126292.2 | 1340935 | EE |
| 3180869 | 2021435 | 534870.4 | 440080.6 | 45283259 | ES |
template("", "maps/recmap1.png", basemap = FALSE, scale = FALSE)
plot.recmap(df, col = NA, border = col, lwd=4, col.text = col)
dev.off()cartog <- recmap(df)
template("", "maps/recmap2.png", basemap = FALSE, scale = FALSE)
plot(cartog[!cartog$name %in% c("IS","MT"),], col = col, border = "white")
dev.off()NB : IS et MT n’ont pu être placés. On les supprime à l’affichage.
cartogramR
library(cartogramR)Le package cartogramR est déveoppé par Pierre-Andre Cornillon et Florent Demoraes de l’université de Rennes. Le package propose les méthodes flow based cartogram (Gastner and Newman 2004), fast flow based cartogram (Gastner, Seguy, and More 2018) et rubber band based cartogram (Dougenik, Chrisman, and Niemeyer 1985). La méthode flow based cartogram est basée sur go_cart.
Attention, cartogramR n’est pas/plus sur le CRAN. L’archive est néanmoins disponible et permet l’installation.
la fonction precartogramR() aide à choisir la taille de la grille de déformation (calcul un peu long). on choisit donc a minima la grille telle que le minimum d’intersections est superieur ou egal a un (ici un pas de grille de 256 peut faire l’affaire)
precarto <- precartogramR(communes, method = "GastnerSeguyMore")
summary(precarto)## L256 L512 L1024 L2048
## Min. 2.00000 5.0000 20.0000 82.000
## 1st Qu. 13.00000 52.0000 208.7500 837.000
## Median 19.00000 77.0000 306.0000 1226.500
## Mean 26.16992 104.6953 418.6816 1674.701
## 3rd Qu. 29.00000 115.2500 460.7500 1850.250
## Max. 405.00000 1613.0000 6469.0000 25884.000
la fonction cartogramR() parmet de calculer le cartogramme en tant que tel selon les 3 méthodes proposées : GastnerSeguyMore" (ou “gsm),”GastnerNewman" (ou “gn),”DougenikChrismanNiemeyer" (ou “dcn”). L’option L=256 permet de choisir la taille de la grille. L’option grid=TRUE (pour les méthodes “gsm” et “gn”) permet d’afficher la grille de déformation. L’option maxit=50 permet de définir le nombre d’itérations max (défaut = 50).
GastnerSeguy <- cartogramR(communes, count="pop2018", method="GastnerSeguyMore", options=list(L=256, grid=TRUE, maxit = 5))## WARNING criterion: 26.455251 > Objective: 0.010000
## Increase maxit or decrease Objective
La fonction make_layer() permet de récupérer la grille de calcul (mais aussi la grille d’origine, les centroides, etc)
grid <- make_layer(GastnerSeguy, type = c("final_graticule"))Ensuite, il est très facile d’afficher le cartogram avec plot (via sf) ou comme précédemment, avec le package mapsf.
template("Gastner, Seguy & More", "maps/gastnerseguy.png", basemap = FALSE, scale = FALSE)
mf_map(GastnerSeguy$cartogram, col = col, border = "white", lwd = 1, add = TRUE)
mf_map(grid, col = NA, border = "#6b4266", lwd = 0.05, add = TRUE)
dev.off()La fonction residuals() permet de calculer les erreurs liées à la déformation (erreur relative : taille finale / taille theorique * 100)
table = cbind(GastnerSeguy$initial_data[,c("id","pop2018")] %>% st_drop_geometry(),
orig_area = GastnerSeguy$orig_area,
final_area = GastnerSeguy$final_area,
errors = residuals(GastnerSeguy, type = "relative error")*100
)
knitr::kable(table[c(0:10),], row.names = F, digits = 1)| id | pop2018 | orig_area | final_area | errors |
|---|---|---|---|---|
| 38001 | 6336 | 27223532 | 39376361.8 | -0.2 |
| 38002 | 1025 | 16353553 | 5636619.1 | -11.7 |
| 38003 | 1127 | 8082831 | 7027719.4 | 0.2 |
| 38004 | 1279 | 10663444 | 7921991.8 | -0.5 |
| 38005 | 954 | 56515131 | 5540617.6 | -6.7 |
| 38006 | 4062 | 34980514 | 24517963.6 | -3.0 |
| 38008 | 28 | 4789504 | 95852.1 | -45.0 |
| 38009 | 1009 | 5121670 | 6220553.5 | -1.0 |
| 38010 | 686 | 13149817 | 4281634.4 | 0.3 |
| 38011 | 1073 | 8596299 | 6774560.2 | 1.4 |
Export au format sf
st_write(as.sf(GastnerSeguy),"files/GastnerSeguy.gpkg", delete_layer=TRUE)
# GastnerSeguy_sf = st_as_sf(table, geometry = GastnerSeguy$cartogram)
# st_write(GastnerSeguy_sf,"files/GastnerSeguy.gpkg", delete_layer=TRUE)Variations
Dots Cartograms
La méthode des dots cartograms est une représentation cartographique à l’intersection des cartogrammes de Dorling et des cartes par points. Les premières cartes utilisant cette méthodes ont été réalisées sur les données du Covid (voir et voir). Article à paraitre. Voir aussi sur Observable
dotcartogram = function(x,var,itermax,onedot,radius){
crs = sf::st_crs(x)
coords <- st_coordinates(st_centroid(sf::st_geometry(x),of_largest_polygon = TRUE))
x <- x[c("id",var)] %>% st_drop_geometry()
x <- data.frame(x, coords)
colnames(x) <- c("id","v","x","y")
x$v <- round(x$v/onedot,0)
x <- x[x$v > 0,]
dots <- x[x$v == 1,c("x","y","v")]
rest <- x[x$v > 1,c("x","y","v")]
nb <- nrow(rest)
for (i in 1:nb){
new <- rest[i,]
new$v <- 1
for (j in 1:rest$v[i]){ dots <- rbind(dots,new)}
}
dots$x <- jitter(dots$x)
dots$y <- jitter(dots$y)
dots$v <- radius
simulation <- circleRepelLayout(x = dots, xysizecols = 1:3,
wrap = FALSE, sizetype = "radius",
maxiter = itermax, weights =1)$layout
circles <- st_buffer(sf::st_as_sf(simulation, coords =c('x', 'y'),
crs = crs), dist = radius)
return(circles)
}onedot = 1000
dc = dotcartogram(x = communes, var = "pop2018", itermax = 120,
onedot = onedot, radius = 600)template(paste0("Un point = ",onedot," personnes"), "maps/dotcartogram.png", scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(dc, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()Anti Cartogramme
Hexbin…
A vous de jouer
Explorez les cartogrammes à l’échelle mondiale.
world <- st_read("data/world_countries.geojson", quiet = TRUE ) %>%
st_transform( "+proj=bertin1953")
world <- world[world$ISO3 != "ATA",]knitr::kable(world[c(0:10),], row.names = F, digits = 1)| ISO2 | ISO3 | ISONUM | NAMEen | NAMEfr | UNRegion | GrowthRate | PopDensity | PopTotal | JamesBond | geometry |
|---|---|---|---|---|---|---|---|---|---|---|
| GQ | GNQ | 226 | Equatorial Guinea | Guinée-Équatoriale | Africa | 2.9 | 30.1 | 845060 | 0 | MULTIPOLYGON (((-514804.2 -… |
| TV | TUV | 798 | Tuvalu | Tuvalu | Oceania | 0.3 | 330.5 | 9916 | 0 | MULTIPOLYGON (((11846963 52… |
| MV | MDV | 462 | Maldives | Maldives | Asia | 1.7 | 1212.2 | 363657 | 0 | MULTIPOLYGON (((5515581 -20… |
| AD | AND | 20 | Andorra | Andorre | Europe | -1.9 | 149.9 | 70473 | 0 | MULTIPOLYGON (((-1013908 15… |
| AG | ATG | 28 | Antigua and Barbuda | Antigua-et-Barbuda | America | 1.0 | 208.7 | 91818 | 0 | MULTIPOLYGON (((-6420562 59… |
| AT | AUT | 40 | Austria | Autriche | Europe | 0.3 | 103.7 | 8544586 | 4 | MULTIPOLYGON (((27563.98 73… |
| BE | BEL | 56 | Belgium | Belgique | Europe | 0.6 | 373.2 | 11299192 | 0 | MULTIPOLYGON (((-621349.5 1… |
| BH | BHR | 48 | Bahrain | Bahren | Asia | 1.4 | 1812.2 | 1377237 | 0 | MULTIPOLYGON (((2804520 -10… |
| BS | BHS | 44 | Bahamas | Bahamas | America | 1.2 | 38.8 | 388019 | 3 | MULTIPOLYGON (((-6563129 25… |
| BB | BRB | 52 | Barbados | Barbade | America | 0.3 | 661.0 | 284215 | 0 | MULTIPOLYGON (((-6530926 70… |
Votre template cartographique
col = "#c291bc"
credits = "Vous, 2021"
theme = mf_theme(x = "default", bg = "#f0f0f0", tab = FALSE,
pos = "center", line = 2, inner = FALSE,
fg = col, mar = c(0,0, 2, 0),cex = 1.9)
templateworld = function(title, file){
mf_export(
world,
export = "png",
width = 1000,
filename = file,
res = 96,
theme = theme,
expandBB = c(-.02,0,-.02,0)
)
mf_title(title)
mf_credits(
txt = credits,
pos = "bottomleft",
col = "#1a2640",
cex = 0.8,
font = 3,
bg = "#ffffff30"
)
}Votre carte de base
templateworld("Le monde en projection Bertin (thx Fil)", "maps/world.png")
mf_map(world, col =col, border = "white", lwd = 0.5, add = TRUE)
dev.off()A vous de jouer…